libs <- c("tidyverse", "sf", "mapSpain", "osmdata", "ggspatial", "dplyr", "readr", "readxl",
"tidyverse", "sf", "ggplot2", "RColorBrewer", "tidygeocoder", "stringi", "leaflet",
"plotly", "dbscan")
installed_libs <- libs %in% rownames(installed.packages())
if (any(installed_libs == F)) {install.packages(libs[!installed_libs])}
invisible(lapply(libs, library, character.only = T))
rm(libs, installed_libs)Tarea 1. Obteniendo datos de OpenStreetMaps con R.
0.Introducción
En este trabajo, analizaremos varios aspectos clave relacionados con el impacto y la distribución de las infraestructuras energéticas y culturales en España, con un enfoque particular en las plantas hidroeléctricas, las zonas industriales en Cataluña y las Fallas en Valencia. A lo largo del estudio, exploraremos cómo estas infraestructuras se distribuyen y cómo afectan al entorno desde diferentes perspectivas geográficas y de agregación.
El análisis comenzará a nivel nacional, donde se estudiarán las plantas hidroeléctricas en España, explorando su localización, la relación con los ríos principales y secundarios, y su impacto en la generación de energía renovable. Posteriormente, se avanzará hacia el nivel de comunidad autónoma (NUTS-2), donde examinaremos las zonas industriales en Cataluña, su relación con la calidad del aire y la contaminación, y cómo la industrialización afecta a las áreas circundantes. Finalmente, en el nivel municipal, nos enfocaremos en las Fallas de Valencia, observando su proximidad a atracciones turísticas y cómo este fenómeno cultural contribuye al paisaje urbano y a la economía local.
Cada una de estas temáticas será explorada en profundidad en los tres niveles de agregación, proporcionando una visión más completa y detallada sobre cómo las infraestructuras energéticas, industriales y culturales impactan en diferentes regiones de España.
Antes de empezar, cargamos los paquetes necesarios a continuación.
Y forzamos que se asuma que hay una conexión a Internet (posiblemente a través de un proxy), independientemente de la configuración real de la red.
curl::has_internet()
assign("has_internet_via_proxy", TRUE, environment(curl::has_internet))
curl::has_internet()
1. Plantas hidroeléctricas en España
Primero, nos centraremos en la extracción y análisis de datos sobre las plantas hidroeléctricas en España, con el objetivo de estudiar su distribución y proximidad a los principales ríos del país. Las plantas hidroeléctricas son una fuente clave de energía renovable en España, y su localización geográfica juega un papel crucial en la eficiencia de la generación de energía.
En esta sección, se pretende no solo analizar la distribución espacial de las plantas hidroeléctricas, sino también explorar su relación con los principales ríos de España, que son los principales canales de la infraestructura hidroeléctrica. Con este análisis, buscaremos obtener información valiosa sobre la capacidad hidroeléctrica del país y sus implicaciones en la planificación energética y la gestión de los recursos hídricos.
1.1 Carga y transformación de los datos
1.1.1 Datos de las plantas hidroeléctricas
Para comenzar, extraemos el polígono que representa la forma de España utilizando el paquete mapSpain. Este objeto, que está en formato sf, se transforma al sistema de coordenadas EPSG:4326, el cual emplea coordenadas geográficas en latitud y longitud. Esta proyección es la misma que utiliza OpenStreetMap (OSM), por lo que mantendremos este sistema de referencia a lo largo de todo el trabajo para asegurar la compatibilidad con los datos obtenidos de OSM.
Una vez transformado el polígono de España, extraemos su bounding box (bbox), que define los límites espaciales del objeto en términos de coordenadas mínimas y máximas. Es importante señalar que OSM no distingue fronteras políticas al extraer datos dentro de un bounding box, por lo que hemos reubicado Canarias dentro del territorio peninsular. Este ajuste evita que el bounding box abarque el norte de África, lo cual podría afectar negativamente al rendimiento del script tanto en términos de almacenamiento como en tiempo de renderizado.
spain <- mapSpain::esp_get_country(moveCAN = T)
spain <- st_transform(spain, 4326)
bbox_spain <- st_bbox(spain) Usando el bbox, cargamos los datos de las plantas hidroeléctricas de España desde OSM con la etiqueta power:plant y plant:source:hydro, de forma que solo nos quedamos con las plantas eléctricas de tipo “hydro”.
power <-osmdata::opq(bbox = bbox_spain) %>%
osmdata::add_osm_feature(key = "power", value = "plant") %>%
osmdata::add_osm_feature(key = "plant:source", value = "hydro") %>%
osmdata::osmdata_sf()Del objeto resultante tipo osmdata_sf extraemos los puntos geométricos a través del atributo osm_points y los guardamos. Comprobamos que su sistema de coordenadas coincide, como esperábamos, con el del bounding box (EPSG:4326).
power_points <-power$osm_points
st_crs(power_points)==st_crs(bbox_spain)[1] TRUEAntes de seguir, echamos un vistazo a los datos: miramos cuántos puntos contiene el objeto y lo visualizamos de forma muy básica.
nrow(power_points)[1] 14483ggplot() +
ggspatial::annotation_map_tile(type = "cartolight") +
geom_sf(data = power_points, color = "black", linewidth = 0.5, alpha = 0.8)
Dado el alto número de puntos según nrow(), en contraste con la cantidad mostrada en el gráfico, es probable que muchos correspondan a vértices de polígonos en lugar de ubicaciones únicas. Esto genera un exceso de puntos duplicados en el mismo lugar. Es crucial eliminarlos para garantizar un recuento preciso, especialmente considerando que los datos han sido recopilados de forma voluntaria.
Primero, replicamos el objeto de los puntos, y lo transformamos al sistema de coordenadas EPSG:3857, que proyecta la superficie del elipsoide terrestre en un plano. Esto permite trabajar con metros como unidad de medida, facilitando un cálculo más preciso de las distancias. Luego, convertimos estas distancias en una matriz, donde la diagonal principal es 0, ya que representa la distancia de cada objeto consigo mismo.
power_points_m <- st_transform(power_points, crs = 3857)
distancias_power <- st_distance(power_points_m)
distancias_power_matrix <- as.matrix(distancias_power)Esta distancia se utiliza como entrada para DBSCAN, un algoritmo de agrupamiento basado en densidad que agrupa puntos según su proximidad. Esto resulta útil en múltiples contextos, y en este caso, nos permite identificar plantas hidroeléctricas sospechosamente cercanas, manteniendo solo un punto representativo de cada polígono.
El parámetro eps define el radio de vecindad, es decir, la distancia máxima en metros dentro de la cual los puntos se consideran parte del mismo grupo. Un valor mayor de eps genera menos grupos. Para este análisis, hemos fijado eps = 1000, asumiendo que las plantas hidroeléctricas diferentes deberían estar separadas al menos por un kilómetro.
El parámetro minPts establece el número mínimo de puntos dentro del vecindario para formar un grupo. En este caso, usamos minPts = 1 para garantizar que los puntos aislados, que no forman parte de un polígono, sean considerados como grupos individuales.
El resultado de DBSCAN asigna un número de grupo a cada punto. Por ejemplo, los puntos octavo, noveno y décimo comparten el grupo 8, lo cual significa que estaban lo suficientemente cerca como para ser considerados parte de la misma estructura. Añadimos la variable que contiene los grupos a los puntos originales, recuperando el sistema de coordenadas EPSG:4326.
dbscan_result_power <- dbscan(distancias_power_matrix, eps = 1000, minPts = 1)
head(dbscan_result_power$cluster,1)[1] 1power_points$cluster <- dbscan_result_power$clusterPara finalizar, solo nos queda agrupar los datos por la variable cluster, eligiendo únicamente un punto por grupo para eliminar el problema de la redundancia. De esta forma hemos reducido el número de puntos casi en la mitad.
power_points_unicos <- power_points %>%
group_by(cluster) %>%
slice(1) %>%
ungroup()
nrow(power_points)[1] 14483nrow(power_points_unicos)[1] 7795Conseguidos los puntos definitivos, decidimos guardar los datos extraídos sobre las zonas industriales en un archivo GeoPackage (GPKG), un formato comúnmente utilizado en sistemas de información geográfica (SIG) para almacenar datos espaciales. GeoPackage permite almacenar diferentes tipos de objetos geoespaciales, como puntos, líneas y polígonos, en un solo archivo, lo que facilita la gestión y el intercambio de datos. A lo largo del trabajo seguiremos realizando este proceso, pero se ha tenido que comentar para evitar el error en el renderizado de “Error: Dataset already exists.” ya que ya existe en nuestro directorio.
#st_write(power_points_unicos, "power.gpkg", layer = "osm_points", driver = "GPKG")
#power_points_unicos <- st_read("power.gpkg", layer = "osm_points")1.1.2 Datos de los ríos
Después de extraer los datos de las plantas hidroeléctricas, es necesario obtener la información sobre los ríos de España. Para ello, utilizamos de nuevo el paquete mapSpain, que nos permite acceder a las geometrías de los ríos en el país. Los ríos están categorizados según su importancia a través del atributo t_rio, donde los ríos principales tienen un valor de 1 y los ríos secundarios un valor de 2.
Dado que los datos abarcan toda la península ibérica, algunos ríos y puntos de plantas hidroeléctricas corresponden a Portugal como vemos en el gráfico anterior. Para evitar incluir información fuera de España, realizamos un filtrado espacial utilizando el polígono del país (spain).
rios_principales <- st_transform(mapSpain::esp_get_rivers() %>% filter(t_rio == 1), 4326)
rios_secundarios <- st_transform(mapSpain::esp_get_rivers() %>% filter(t_rio == 2), 4326)
rios_principales_filtrados <- st_filter(rios_principales, spain)
rios_secundarios_filtrados <- st_filter(rios_secundarios, spain)
power_unique_filtrados <- st_filter(power_points_unicos, spain)Para saber si las plantas hidroeléctricas se encuentran en un río principal o secundario, filtramos aquellas que se encuentran a una distancia máxima de 500 metros de estos ríos. Para ello, utilizamos la función st_is_within_distance(), que evalúa qué puntos de power_unique_filtrados están dentro de la distancia especificada con respecto a las líneas que representan los ríos.
Esta función devuelve una lista donde cada elemento representa una planta hidroeléctrica y contiene una referencia a los ríos cercanos dentro de la distancia especificada (500 metros). Con esta lista identifica y filtramos aquellas plantas que están cerca de ríos. Ya que lengths() calcula la cantidad de ríos cercanos a cada planta, buscamos aquellos con un valor mayor a cero.
power_rios_principales <- st_is_within_distance(power_unique_filtrados, rios_principales_filtrados, dist = 500)
power_rios_principales <- power_unique_filtrados[lengths(power_rios_principales) > 0, ]
power_rios_secundarios <- st_is_within_distance(power_unique_filtrados, rios_secundarios_filtrados, dist = 500)
power_rios_secundarios <- power_unique_filtrados[lengths(power_rios_secundarios) > 0, ]1.2 Graficado e interpretación de los resultados
Finalmente, utilizamos ggplot() para visualizar los resultados. Sobre un mapa base con estilo “cartolight” se representan cuatro capas: primero, los ríos secundarios con un color azul claro y una línea más delgada y transparente (linewidth = 0.5, alpha = 0.8) para diferenciarlos visualmente de los ríos principales, que se trazan con el mismo color pero con una línea más gruesa (linewidth = 1). Y segundo, las plantas hidroeléctricas cercanas a ríos secundarios, que se representan como puntos azules de menor tamaño (size = 0.5), mientras que las cercanas a ríos principales se destacan con un tamaño mayor (size = 1).
ggplot() +
ggspatial::annotation_map_tile(type = "cartolight", zoom = 8) +
geom_sf(data = rios_secundarios_filtrados, color = "lightblue", linewidth = 0.5, alpha = 0.8) +
geom_sf(data = rios_principales_filtrados, color = "lightblue", linewidth = 1) +
geom_sf(data = power_rios_secundarios, color = "blue", size = 0.5) +
geom_sf(data = power_rios_principales, color = "blue", size = 1) +
theme_minimal() +
labs(title = "Ubicación de las plantas hidroeléctricas sobre los ríos principales")
Como podemos observar, los ríos con más plantas hidroeléctricas parecen ser el Tajo y la zona más norteña del Ebro. Esto sugiere que estos tienen una mayor capacidad para la generación de energía hidroeléctrica, probablemente debido a su caudal, longitud y desnivel en ciertos tramos.
El río Tajo, el más largo de la península ibérica, atraviesa varias regiones de España y cuenta con múltiples embalses y presas a lo largo de su recorrido, lo que facilita la instalación de infraestructuras hidroeléctricas. Por otro lado, el Ebro, especialmente en su zona norte, recorre una orografía más montañosa, lo que favorece el aprovechamiento de la energía potencial del agua para la generación de electricidad.
Este patrón también podría estar influenciado por factores históricos y políticos, como la planificación energética en España, que ha priorizado determinadas cuencas hidrográficas para el desarrollo de la energía hidroeléctrica.
1.3 Comparación datos de fuente alternativa
Dando por hecho que las grandes presas contienen las centrales hidroeléctricas, y estas se encuentran en grandes ríos (principales y secundarios), contamos el número de plantas que se encuentran en rios principales y las sumamos a las que hay en los secundarios.
nrow(power_rios_principales) + nrow(power_rios_secundarios)[1] 1764Según el Ministerio para la Transición Ecológica y el Reto Demográfico, “En la actualidad el número de grandes presas supera las 1.200 con una capacidad aproximada de 56.000 hm3.” Por lo tanto, el número de plantas hidroeléctricas que encontramos en los ríos principales y secundarios (1761) es un valor cercano a la realidad en términos de la infraestructura hidroeléctrica en España.
De todas formas, la discrepancia entre el número de plantas hidroeléctricas reportadas en OSM y las cifras del Ministerio puede deberse a varios factores. OSM incluye todas las plantas hidroeléctricas, sin distinguir entre grandes presas y pequeñas instalaciones, lo que podría dar como resultado un número superior de plantas, incluyendo aquellas de menor capacidad que no están reflejadas en el Ministerio. Además, OSM se nutre de datos proporcionados por voluntarios, lo que implica que puede contener información más detallada o amplia, especialmente en áreas menos cubiertas por datos oficiales, o por el contrario, información menos fiable.
El hecho de que tengamos un total de más de 1.200 plantas en estos ríos refleja la importancia de los ríos principales y secundarios en la generación de energía hidroeléctrica en el país. Además, la capacidad total de almacenamiento de agua (56.000 hm³) demuestra la magnitud del potencial hidroeléctrico del país. Este volumen de agua almacenada no solo permite la generación de energía en períodos de alta demanda, sino que también juega un papel clave en la gestión del recurso hídrico en España, especialmente en tiempos de sequía.
2. Industrias en Cataluña
A continuación, vamos a analizar cómo las zonas industriales en Cataluña impactan en la calidad del aire de las áreas circundantes. La industrialización, si bien es un factor clave para el crecimiento económico y el desarrollo, genera externalidades negativas como la contaminación del aire, que puede afectar tanto a la salud de las personas como al entorno natural. En este contexto, las emisiones de gases contaminantes provenientes de fábricas, centrales eléctricas y otras instalaciones industriales se dispersan en el aire y pueden afectar a poblaciones cercanas, especialmente en zonas urbanas densamente pobladas.
El análisis consistirá en estudiar la distribución geográfica de las zonas industriales en Cataluña y la calidad del aire en estas áreas. Este análisis permitirá identificar las áreas más afectadas por la contaminación industrial y su relación con la proximidad a los centros de actividad industrial.
El objetivo final es aportar información sobre los efectos de la industrialización sobre la calidad del aire, y proporcionar información que permita a los responsables políticos tomar decisiones informadas para mitigar estos impactos negativos, equilibrando el desarrollo económico con la protección del medio ambiente y la salud pública.
2.1 Carga y transformación de datos
2.1.1 Datos de calidad del aire
Por un lado, hemos descargado del Ministerio para la Transición Ecológica y el Reto Demográfico un archivo titulado “Datos oficiales Calidad del Aire 2023”. Este archivo contiene, entre otras cosas, información detallada sobre la concentración de grandes partículas presentes en el aire en diversas estaciones de monitoreo situadas a lo largo de Cataluña. Las estaciones de calidad del aire registran continuamente datos sobre contaminantes atmosféricos, tales como partículas en suspensión (PM10 y PM2.5), que son especialmente relevantes para evaluar el impacto de la contaminación industrial y de tráfico en la salud pública. Con estos datos, podemos analizar cómo varían los niveles de contaminantes en las zonas cercanas a las industrias, lo que nos ayudará a entender mejor la relación entre la localización de las actividades industriales y la calidad del aire en la región.
PM10 <- read_delim("PM10_DD_2023.csv",
delim = ";", escape_double = FALSE, trim_ws = TRUE)
head(PM10)# A tibble: 6 × 38
PROVINCIA MUNICIPIO ESTACION MAGNITUD PUNTO_MUESTREO ANNO MES D01 D02
<dbl> <dbl> <dbl> <dbl> <chr> <dbl> <dbl> <dbl> <dbl>
1 1 59 9 10 01059009_10_M 2023 1 12 6
2 1 59 9 10 01059009_10_M 2023 2 5 13
3 1 59 9 10 01059009_10_M 2023 3 NA NA
4 1 59 9 10 01059009_10_M 2023 4 NA NA
5 1 59 9 10 01059009_10_M 2023 5 7 10
6 1 59 9 10 01059009_10_M 2023 6 16 14
# ℹ 29 more variables: D03 <dbl>, D04 <dbl>, D05 <dbl>, D06 <dbl>, D07 <dbl>,
# D08 <dbl>, D09 <dbl>, D10 <dbl>, D11 <dbl>, D12 <dbl>, D13 <dbl>,
# D14 <dbl>, D15 <dbl>, D16 <dbl>, D17 <dbl>, D18 <dbl>, D19 <dbl>,
# D20 <dbl>, D21 <dbl>, D22 <dbl>, D23 <dbl>, D24 <dbl>, D25 <dbl>,
# D26 <dbl>, D27 <dbl>, D28 <dbl>, D29 <dbl>, D30 <dbl>, D31 <dbl>Este conjunto de datos ha sido agrupado por PROVINCIA, MUNICIPIO y ESTACION, calculando la media de la variable D01, que representa uno de los tipos de partículas contaminantes en el aire. De esta forma, hemos obtenido un valor representativo de la calidad del aire en cada localidad o zona donde se encuentran las estaciones de monitoreo. Además, hemos cargado los metadatos de las estaciones, los cuales contienen las coordenadas exactas de cada una. Esta información es fundamental, ya que nos permitirá ubicar las estaciones en el mapa y relacionar visualmente los niveles de contaminación con las áreas geográficas.
PM10_resumido <- PM10 %>%
group_by(PROVINCIA, MUNICIPIO, ESTACION) %>% # Sin comillas
summarise(D01 = mean(D01, na.rm = TRUE), .groups = "drop") # Resumir valores
Meta <- read_excel("Metainformacion2023_19022025.xlsx", sheet="Estaciones2023")
head(Meta)# A tibble: 6 × 21
COD_LOCAL PROVINCIA MUNICIPIO ESTACION COD_ESTACION_DEM N_RED NOMBRE
<dbl> <dbl> <dbl> <dbl> <chr> <chr> <chr>
1 1022001 1 22 1 ES1672A CCAA País Vasco EL CI…
2 1036004 1 36 4 ES1349A CCAA País Vasco LLODIO
3 1051001 1 51 1 ES1544A CCAA País Vasco AGURA…
4 1055001 1 55 1 ES1489A CCAA País Vasco VALDE…
5 1059008 1 59 8 ES1502A CCAA País Vasco AVENI…
6 1059009 1 59 9 ES1492A CCAA País Vasco TRES …
# ℹ 14 more variables: FECHA_INI <dttm>, FECHA_FIN <dttm>, LATITUD_G <dbl>,
# LONGITUD_G <dbl>, ALTITUD <dbl>, N_CCAA <chr>, N_PROVINCIA <chr>,
# N_MUNICIPIO <chr>, TIPO_ESTACION <chr>, TIPO_AREA <chr>,
# TIPO_SUBAREA_RURAL <chr>, EST <chr>, ZONA <chr>, DIRECCION <chr>Después, hemos realizado una unión de ambas bases de datos utilizando un left_join para asegurar que todas las estaciones con información sobre las partículas (calculada previamente) queden asociadas a los metadatos correspondientes. Esto se ha hecho a través de una clave combinada única que se genera al unir los campos PROVINCIA, MUNICIPIO y ESTACION.
datos_completos <- left_join(PM10_resumido, Meta, by = c("PROVINCIA", "MUNICIPIO", "ESTACION"))
head(datos_completos)# A tibble: 6 × 22
PROVINCIA MUNICIPIO ESTACION D01 COD_LOCAL COD_ESTACION_DEM N_RED NOMBRE
<dbl> <dbl> <dbl> <dbl> <dbl> <chr> <chr> <chr>
1 1 59 9 9.7 1059009 ES1492A CCAA Paí… TRES …
2 3 2 2 22.2 3002002 ES1847A CCAA Com… AGOST
3 3 9 6 12.2 3009006 ES1623A CCAA Com… ALCOI…
4 3 14 6 22.8 3014006 ES1635A CCAA Com… ALACA…
5 3 65 6 19.2 3065006 ES1624A CCAA Com… ELX-A…
6 3 65 7 18.8 3065007 ES1849A CCAA Com… ELX-P…
# ℹ 14 more variables: FECHA_INI <dttm>, FECHA_FIN <dttm>, LATITUD_G <dbl>,
# LONGITUD_G <dbl>, ALTITUD <dbl>, N_CCAA <chr>, N_PROVINCIA <chr>,
# N_MUNICIPIO <chr>, TIPO_ESTACION <chr>, TIPO_AREA <chr>,
# TIPO_SUBAREA_RURAL <chr>, EST <chr>, ZONA <chr>, DIRECCION <chr>Una vez filtrados los datos para centrarnos únicamente en las estaciones de la comunidad autónoma de Cataluña, creamos un objeto espacial sf (simple features) con las coordenadas geográficas de las estaciones. Para ello, utilizamos las columnas LONGITUD_G y LATITUD_G, que contienen las ubicaciones exactas de cada estación de monitoreo de calidad del aire.
datos <- datos_completos %>% filter(N_CCAA == "CATALUÑA") # Filtramos solo Barcelona
estaciones_sf <- st_as_sf(datos, coords = c("LONGITUD_G", "LATITUD_G"), crs = 4326)
head(estaciones_sf)Simple feature collection with 6 features and 20 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 2.1331 ymin: 41.3788 xmax: 2.2507 ymax: 41.7711
Geodetic CRS: WGS 84
# A tibble: 6 × 21
PROVINCIA MUNICIPIO ESTACION D01 COD_LOCAL COD_ESTACION_DEM N_RED NOMBRE
<dbl> <dbl> <dbl> <dbl> <dbl> <chr> <chr> <chr>
1 8 14 1 20.0 8014001 ES2088A CCAA Cat… AIGUA…
2 8 15 23 19.5 8015023 ES2027A CCAA Cat… BADAL…
3 8 19 4 19.8 8019004 ES0691A CCAA Cat… BARCE…
4 8 19 28 24.1 8019028 ES0559A CCAA Cat… BARCE…
5 8 19 42 22.3 8019042 ES1396A CCAA Cat… BARCE…
6 8 19 43 22.1 8019043 ES1438A CCAA Cat… BARCE…
# ℹ 13 more variables: FECHA_INI <dttm>, FECHA_FIN <dttm>, ALTITUD <dbl>,
# N_CCAA <chr>, N_PROVINCIA <chr>, N_MUNICIPIO <chr>, TIPO_ESTACION <chr>,
# TIPO_AREA <chr>, TIPO_SUBAREA_RURAL <chr>, EST <chr>, ZONA <chr>,
# DIRECCION <chr>, geometry <POINT [°]>2.1.2 Datos de las industrias
Por otro lado, utilizamos OSM para obtener la información geométrica de Cataluña, lo que nos permitirá identificar las zonas industriales dentro de la comunidad autónoma. A través de una consulta a la base de datos de OSM, extraemos objetos de tipo industrial utilizando la etiqueta Key:industrial, la cual según OSMWiki, hace referencia a “un tipo de industria o un tipo de objeto industrial”. Esta etiqueta nos proporciona información detallada sobre las zonas industriales en Cataluña, como fábricas, almacenes y otros complejos industriales.
bbox_cat <- osmdata::getbb("Cataluña")
q <- bbox_cat %>%
osmdata::opq() %>%
osmdata::add_osm_feature(key = "industrial")
ind <- osmdata::osmdata_sf(q)Extraemos las capas correspondientes a puntos y líneas desde el objeto obtenido de OSM, almacenándolas en las variables ind_points y ind_lines, respectivamente. Luego, utilizamos la función st_write() para guardar estas capas en un archivo GeoPackage llamado “ind.gpkg”, especificando las capas y el formato correspondiente.
ind_points <- ind$osm_points
ind_lines <- ind$osm_lines
#st_write(ind_points, "ind.gpkg", layer = "osm_points", driver = "GPKG")
#st_write(ind_lines, "ind.gpkg", layer = "osm_lines", driver = "GPKG")
#ind_points <- st_read("ind.gpkg", layer = "osm_points")
#ind_lines <- st_read("ind.gpkg", layer = "osm_lines")2.2 Graficado e interpretación de los resultados
Con todos los datos necesarios cargados, podemos proceder a la visualización. Para ello, creamos los buffers alrededor de cada estación de calidad del aire, los cuales son círculos cuyo centroide serán las estaciones. Primero, transformamos las coordenadas de las estaciones al sistema de coordenadas EPSG:3857, que es proyectado y usa metros como unidades. A continuación, calculamos el tamaño de los buffers, multiplicando los valores de D01 por 1000 metros. Este tamaño está relacionado con la concentración de partículas en el aire en cada estación, lo que nos da una idea del área afectada por la contaminación. Generamos los buffers utilizando la función st_buffer() de sf, donde cada estación se convierte en un área circular cuyo radio depende del valor de D01.
estaciones_sf_m <- st_transform(estaciones_sf, crs = 3857)
buffer_size <- estaciones_sf_m$D01* 1000
estaciones_sf$buffer_size <- buffer_size
estaciones_buffer <- estaciones_sf %>%
st_buffer(dist = estaciones_sf$buffer_size) # Crea los buffers usando el tamaño definidoEste gráfico nos muestra cinco capas. Para empezar, por debajo queda el mapa de Cataluña con el tipo o estilo “osm”. Después le hemos añadido una capa de los buffers de las estaciones, los cuales hemos pintado según su variacartolightble D01 con la escala “magma” invertida, de forma que las zonas con aire más contaminado se vean más oscuras, además, esta variable también se ve con el tamaño de cada buffer, que da más visibilidad a las zonas contaminadas. Por último, hemos añadido tres capas que representan los tres niveles de la geometría que hemos extraido de los “objetos industriales”: puntos, lineas y polígonos, que hemos pintado de color negro para que resalten sobre los buffers.
Cabe destacar que en este caso no nos hemos preocupado por la redundancia de datos ya que solo nos preocupa su graficado, que se vería poco afectado.
ggplot() +
ggspatial::annotation_map_tile(type = "osm", zoom = 10) +
geom_sf(data = estaciones_buffer, aes(fill = D01), color = NA, alpha = 0.3) +
geom_sf(data = ind_points, color = "gray7", size=1) +
geom_sf(data = ind_lines, color = "black", size=1) +
scale_fill_viridis_c(name = "D01", option = "magma", direction = -1) +
scale_color_identity() +
theme_void() +
labs(title = "Relación entre partículas D01 y situación de las industrias catalanas")
Al observar el gráfico, podemos notar que hay una aparente relación entre las zonas más industrializadas (habitualmente las grandes ciudades) y las zonas con peor calidad del aire, ya que las áreas más oscuras (de mayor concentración de partículas) corresponden a estas Sin embargo, también es evidente que la calidad del aire en áreas rurales no está tan bien representada, lo que nos deja con la duda de si la falta de industria en esas zonas realmente conlleva una mejor calidad del aire. Además, en ciudades como Barcelona, la visualización se ve saturada por la cantidad de información, lo que dificulta un análisis claro de la situación.
Para una investigación más precisa, sería necesario contar con más datos sobre la calidad del aire en zonas rurales y un sistema de zonificación que permita estudiar el nivel de contaminación de manera más segmentada. De este modo, sería posible obtener conclusiones más fiables sobre la relación entre la industrialización y la calidad del aire.
3. Fallas en Valencia
A continuación, vamos a analizar la relación entre las Fallas y las principales atracciones turísticas de Valencia. Las Fallas, uno de los eventos más importantes de la ciudad, atraen a miles de turistas cada año, lo que genera un impacto significativo tanto en la economía local como en la vida cultural de la ciudad.
Este análisis consistirá en estudiar la proximidad de las Fallas a los museos y monumentos más relevantes de la ciudad, evaluando cómo las Fallas pueden contribuir positivamente a la generación de turismo cultural. Utilizando herramientas geoespaciales, se identificará qué Fallas están más cerca de estos puntos turísticos y se analizará la distribución de estos espacios a lo largo de la ciudad.
El objetivo es proporcionar una visión clara sobre la relación entre las Fallas y los principales atractivos turísticos, con el fin de ayudar a la planificación urbana y turística en Valencia, optimizando los flujos de visitantes y garantizando que tanto el evento como las atracciones culturales puedan ser disfrutados de manera eficiente y sostenible.
3.1 Carga y transformación de los datos
3.1.1 Datos de las fallas
Primero hemos creado un dataframe llamado fallas_data que contiene información sobre las fallas de Valencia. Específicamente, tenemos dos columnas: una con el nombre de cada falla y otra con su localización exacta en formato de dirección. Estas direcciones las normalizamos convirtiendo todas las letras a mayúsculas y quitando los caracteres especiales (como tildes) con la función stri_trans_general. Esto es útil porque el geocodificador puede no funcionar correctamente si la dirección tiene caracteres especiales o tildes.
fallas_data <- data.frame(
Falla = c(
"Falla L'antiga de Campanar",
"Exposició-Misser Mascó",
"Falla Plaza del Pilar",
"Falla Convento Jerusalén",
"Falla Na Jordana",
"Falla Sueca Literato Azorín",
"Falla Reino de Valencia - Duc de Calabria",
"Falla Cuba Literato Azorín",
"Falla Almirante Cadarso - Conde Altea",
"Falla Espartero",
"Mestre Gozalbo-Comte d'Altea",
"Falla Grabador Esteve - Cirillo Amoros",
"Falla Arzobispo Olaechea - San Marcelino",
"Falla Islas Canarias Trafalgar Samuel Ros",
"Falla Linterna Na Robella",
"Falla Santa Genoveva Torres - Arquitecto Tolsa-Alfahuir - La Nova d'Orriols",
"Falla Quart i Palomar",
"Falla Ciscar Borriana",
"Falla Jeronima Gales - Litografo Pascual i Abad",
"Falla San Vicente - Periodista Azzati",
"Falla Comte de Salvatierra-Ciril Amorós",
"Falla Just Vilar-Mercat del Cabanyal",
"Falla Ribera i Convent de Santa Clara",
"Falla Mercado Central",
"Falla Quart Extramurs",
"Falla Plaza de la Merced"
),
Localizacion = c(
"Carrer d'Aparicio Albiñana, 2, 46015 Valencia",
"Carrer del Naturalista Arévalo Baca, 46010 València, Valencia",
"Plaça del Pilar, 46001 València",
"Calle Convento Jerusalén, 26, 46007 Valencia",
"Carrer de Salvador Giner, 9, 46003 València, Valencia",
"Carrer de Sueca, 50, 46004 València, Valencia",
"C. del Duque de Calabria, 6, 46005 Valencia",
"Calle Literato Azorín, 39, 46004 València",
"Calle del Almirante Cadarso, 19, 46005 Valencia",
"C/ del Pare Jofré, 24",
"Calle del Comte d'Altea, 5",
"Carrer del Gravador Esteve, València",
"C/ Arzobispo Olaechea, 30, 46017 València",
"Carrer de les Illes Canàries, 108, 46023 València",
"C/ de la Beata, 8, Ciutat Vella, 46001 València, Valencia",
"C. de Santa Genoveva Torres, 30, 46019 València",
"C. Quart, 29, 46001 Valencia",
"Carrer de Ciscar, 36, 46005 València",
"C. del Músic Penella, 10, 46017 València",
"Carrer del Periodista Azzati, 46002 València",
"C/ del Cronista Carreres, 7, 46003",
"Carrer de Francesc Baldomar, 84, 46011 València",
"Carrer de Ribera, 8, 46002 València",
"Plaça del Mercat, 46001 València, Valencia",
"Carrer de la Democràcia, 74, 46018 València",
"Plaça de la Mercè, 46001 València, Valencia"
)
)
fallas_data$Localizacion <- toupper(stri_trans_general(fallas_data$Localizacion,"Latin-ASCII"))
head(fallas_data) Falla
1 Falla L'antiga de Campanar
2 Exposició-Misser Mascó
3 Falla Plaza del Pilar
4 Falla Convento Jerusalén
5 Falla Na Jordana
6 Falla Sueca Literato Azorín
Localizacion
1 CARRER D'APARICIO ALBINANA, 2, 46015 VALENCIA
2 CARRER DEL NATURALISTA AREVALO BACA, 46010 VALENCIA, VALENCIA
3 PLACA DEL PILAR, 46001 VALENCIA
4 CALLE CONVENTO JERUSALEN, 26, 46007 VALENCIA
5 CARRER DE SALVADOR GINER, 9, 46003 VALENCIA, VALENCIA
6 CARRER DE SUECA, 50, 46004 VALENCIA, VALENCIASobre este dataframe usamos la función geo de tidygeocoder para realizar la geocodificación de las direcciones. Esto significa que convertimos las direcciones en coordenadas geográficas (latitud y longitud). En este caso, no especificamos el proveedor del servicio, de forma que utiliza por defecto el método “osm” o Nominatim, un servicio de geocodificación gratuito basado en datos de OpenStreetMap.
fallas_crds <- geo(street = fallas_data$Localizacion,
country = rep("Spain",nrow(fallas_data)))
sum(is.na(fallas_crds$lat))[1] 19Como vemos, el problema de este servicio es que no suele ser muy efectivo, a pesar de que le demos información adicional, como en este caso que hemos añadido country = “Spain” para que no fuera a ningún otro país con direcciones similares.
Utilizamos por lo tanto, un método alternativo de geocodificación: el servicio de ArcGIS, que ofrece un servicio alternativo más preciso y que puede manejar mejor estas ambigüedades. Este servicio es conocido por ser más robusto para algunas localizaciones y direcciones complicadas como es el caso.
fallas_crds_arcgis <- tidygeocoder::geocode(fallas_data, address = Localizacion, method = "arcgis")
fallas_sf <- sf::st_as_sf(fallas_crds_arcgis, coords=c("long","lat"), crs=4326, na.fail=F)
sum(is.na(fallas_crds_arcgis$lat))[1] 1De esta forma, encontramos que solo una de las fallas (Falla Espartero) no se geocodificó correctamente. Entonces, ajustamos manualmente la ubicación utilizando las coordenadas correctas. Específicamente, asignamos un punto con las coordenadas de latitud 39.47 y longitud -0.38 que extraemos de una simple búsqueda en GoogleMaps. Lo hacemos de forma manual porque un único valor no conlleva a penas esfuerzo.
Cabe mencionar que a la hora de convertir el objeto en un objeto sf, el orden de longitud y latitud es esencial. Por ejemplo, en Alboraya, la latitud es 39.5º y la longitud -0.35, mientras que en Irangi, Kenia, la longitud es 39.5 y la latitud -0.35.
fallas_sf$geometry[10]Geometry set for 1 feature (with 1 geometry empty)
Geometry type: POINT
Dimension: XY
Bounding box: xmin: NA ymin: NA xmax: NA ymax: NA
Geodetic CRS: WGS 84fallas_sf$geometry[10] <- st_sfc(st_point(c(-0.38, 39.47)), crs = 4326)
fallas_sf$geometry[10]Geometry set for 1 feature
Geometry type: POINT
Dimension: XY
Bounding box: xmin: -0.38 ymin: 39.47 xmax: -0.38 ymax: 39.47
Geodetic CRS: WGS 843.1.2 Datos de las atracciones turísticas principales
A continuación, trabajamos con un conjunto de datos que contiene información sobre las principales atracciones turísticas en Valencia, como museos y monumentos.
Para ello, primero definimos un área de interés en Valencia mediante un bounding box, el cual hemos especificado manualmente. Usando este área de interés, consultamos los datos de OpenStreetMap para obtener información sobre las atracciones turísticas de la ciudad, centrándonos específicamente en museos y monumentos. A partir de los datos obtenidos, extraemos los puntos y polígonos correspondientes a las atracciones turísticas. Además, de los polígonos, calculamos los centroides, los cuales también representan puntos geográficos.
bbox_valencia <- c(xmin = -0.41, ymin = 39.43, xmax = -0.25, ymax = 39.51)
bbox_valencia <- sf::st_bbox(bbox_valencia, crs = 4326)
turism <-osmdata::opq(bbox = bbox_valencia) %>%
osmdata::add_osm_feature(key = "tourism", value = c("museum", "monument")) %>%
osmdata::osmdata_sf()
turism_points <- turism$osm_points
turism_pol <- turism$osm_polygons
turism_centroids <- st_centroid(turism_pol)Primero, en cuanto a los puntos, resolvemos el problema de la redundancia de nuevo. Calculamos las distancias entre los puntos de atracciones turísticas utilizando la función st_distance, tras haber duplicado el objeto pero con el sistema de coordenadas EPSG:3857 que usa las unidades de medida de megtros.
Luego, identificamos los puntos “solos”, es decir, aquellos que están a más de 100 metros de cualquier otro punto de atracción turística. Para lograrlo, utilizamos un umbral de 100 metros. Primero, calculamos la matriz de distancias entre todos los puntos de atracciones turísticas. Después, para cada punto, verificamos si existe otro punto a menos de 100 metros de distancia. Esto se hace al revisar todas las distancias en la fila correspondiente a cada punto (sin contar la distancia del punto consigo mismo, que está en la diagonal de la matriz de distancias). Si todas las distancias a otros puntos son mayores a 100 metros, consideramos que ese punto está “solo”.
Finalmente, filtramos los puntos que cumplen con esta condición y los almacenamos en turism_points_solos. Estos puntos representan las atracciones turísticas que están aisladas de las demás en cuanto a distancia.
turism_points_m <- st_transform(turism_points, crs = 3857)
distancias <- st_distance(turism_points_m)
solos <- apply(distancias, 1, function(x) all(x[-which.min(x)] > 100))
turism_points_solos <- turism_points[solos, ]Segundo, en cuanto a los centroides, nos encontramos con el probolema de que algunos edificios se componen por diversos polígonos (dos edificios de un mismo museo, patios interiores, etc). Para solucionarlo, volvemos a realizar un análisis de agrupación utilizando el algoritmo DBSCAN.
Primero, volvemos a transformar las coordenadas de los centroides de los polígonos de atracciones turísticas al sistema de coordenadas proyectado EPSG: 3857 para realizar cálculos de distancia más precisos en unidades de metros.
Luego, calculamos la distancia entre los centroides utilizando la función st_distance, que nos devuelve una matriz de distancias entre todos los puntos de los centroides, y convertimos esta matriz.
Con la matriz de distancias, aplicamos el algoritmo DBSCAN para agrupar los centroides. El parámetro eps define el radio de búsqueda para los puntos cercanos (en este caso, 1000 metros) y minPts especifica el número mínimo de puntos requeridos para formar un grupo (en este caso, 1, es decir, cada punto es considerado como un grupo por sí mismo, si no encuentra otros puntos cercanos). El resultado de DBSCAN asigna un número de cluster a cada centroide.
Posteriormente, seleccionamos solo un punto por grupo con slice(1). Este paso nos permite obtener un punto representativo para cada grupo de centroides cercanos.
turism_centroids_m <- st_transform(turism_centroids, crs = 3857)
distancias_centroids <- st_distance(turism_centroids_m)
distancias_centroids_matrix <- as.matrix(distancias_centroids)
dbscan_result <- dbscan(distancias_centroids_matrix, eps = 1000, minPts = 1)
turism_centroids$cluster <- dbscan_result$cluster
centroids_unicos <- turism_centroids %>%
group_by(cluster) %>%
slice(1) %>%
ungroup()Finalmente, combinamos estos puntos representativos con aquellos puntos “solos” anteriores, creando un conjunto final de puntos únicos que representan las atracciones turísticas. Hemos seleccionado previamente las columnas geometry y name para que sus dimensiones coincidieran y pudiéramos realizar el rbind.
centroids_unicos <- centroids_unicos %>%
select(geometry, name)
turism_points_solos <- turism_points_solos %>%
select(geometry, name)
todos_museos <- rbind(turism_points_solos, centroids_unicos)
head(todos_museos)Simple feature collection with 6 features and 1 field
Geometry type: POINT
Dimension: XY
Bounding box: xmin: -0.3849821 ymin: 39.46628 xmax: -0.3746678 ymax: 39.48716
Geodetic CRS: WGS 84
name
1981174882 <NA>
2508173010 Museo Taurino
2700282752 Museo de los Soldaditos de Plomo L'Iber
5100017480 Bombas Gens Centre d'Art
5126811958 Museo Nacional de Cerámica y Artes Suntuarias González Martí
5596900823 <NA>
geometry
1981174882 POINT (-0.3806511 39.48716)
2508173010 POINT (-0.3756447 39.46628)
2700282752 POINT (-0.3782275 39.47672)
5100017480 POINT (-0.3849821 39.48555)
5126811958 POINT (-0.3746678 39.47252)
5596900823 POINT (-0.3760315 39.47458)Para una base de datos más completa, hemos añadido manualmente aquellos nombres de los museos que no se registraron en su momento correctamente simplemente buscando en GoogleMaps qué museos se encontraban en cada coordenada correspondiente.
todos_museos[1,1] <- "Exposició: Espai Ferroviari a Marxalenes"
todos_museos[6,1] <- "Set Espai d'Art"
todos_museos[13,1] <- "Museu de Belles Arts de València"
todos_museos[15,1] <- "Museo de las Ciencias"
todos_museos[17,1] <- "Casa Museo Benlliure"
todos_museos[19,1] <- "Museo Palacio del Marqués de Campo"
todos_museos[22,1] <- "Museo del Corpus - Casa de las Rocas"
todos_museos[25,1] <- "Museo Histórico Militar"De nuevo guardamos el objeto en el directorio con la extensión .gpkg para asegurar su reproducibilidad.
#st_write(todos_museos, "museos.gpkg", layer = "osm_points", driver = "GPKG")
#power_points_unicos <- st_read("power.gpkg", layer = "osm_points")3.2 Graficado e interpretación de los resultados
Para graficar, creamos de nuevo un buffer alrededor de cada falla. Usamos la función st_intersects() para verificar si los buffers de las fallas (las zonas de influencia) intersectan con los puntos que representan los museos. Esto devuelve una lista de índices que indica qué museos están dentro de cada buffer de falla.
Utilizamos la función sapply() para calcular la cantidad de museos dentro de cada buffer de falla. Este número se guarda en una nueva columna museos_buffer dentro del dataframe fallas_sf.
todos_museos_m <- st_transform(todos_museos, crs = 3857)
fallas_sf_m <- st_transform(fallas_sf, crs = 3857)
fallas_buffer <- st_buffer(fallas_sf_m, dist = 700)
museos_en_buffer <- st_intersects(fallas_buffer, todos_museos_m)
fallas_sf$museos_buffer <- sapply(museos_en_buffer, length)
fallas_sf$museos_buffer [1] 0 1 2 2 4 0 0 0 0 3 1 0 0 0 5 0 2 0 0 5 2 0 5 6 1 8Esta vez generamos un mapa interactivo utilizando leaflet ya que tiene más interés.
Añadimos un mapa base (fondo) de OpenStreetMap utilizando la capa CartoDB.Positron; agregamos los museos al mapa como puntos (círculos de color naranja), con etiquetas que muestran el nombre de cada museo; agregamos las fallas al mapa cuyo color está determinado por el número de museos cercanos (usando la paleta de colores que creamos antes con pal), además, se incluyen ventanas emergentes (popups) que muestran el nombre de la falla y el número de atracciones turísticas cercanas.
paleta_colores <- colorNumeric(palette = c("yellow", "darkgreen"),
domain = fallas_sf$museo_count)
pal <- colorNumeric(palette = "viridis", domain = fallas_sf$museos_buffer, reverse =TRUE)
leaflet() %>%
addProviderTiles("CartoDB.Positron") %>%
addCircleMarkers(data = todos_museos,
color = "orange",
radius = 8,
stroke = FALSE,
fillOpacity = 0.7,
label = ~name,
group = "Museos/Monumentos") %>%
addCircleMarkers(data = fallas_sf,
color = ~pal(museos_buffer),
radius = 10,
stroke = FALSE,
fillOpacity = 0.7,
popup = ~paste(Falla, "<br>",
"Atracciones turísticas cercanas: ", museos_buffer),
group = "Fallas") %>%
addLegend(position = "topright",
pal = pal,
values = fallas_sf$museos_buffer, # Valores de museos cercanos
title = paste("Número de atracciones <br> turísticas cercanas"),
labFormat = labelFormat(digits = 0), # Sin decimales
opacity = 1)El mapa es especialmente interesante porque muestra no solo la ubicación de cada elemento, sino también las interacciones entre ellos: cómo las fallas están influenciadas o relacionadas con la cantidad de museos cercanos. Esta representación ayuda a identificar puntos de concentración turística, resaltar áreas estratégicas para eventos culturales, e incluso planificar futuros proyectos turísticos o urbanos. Además, las ventanas emergentes añaden una capa de interactividad, permitiendo explorar datos detallados de forma atractiva y dinámica. Es una herramienta poderosa para entender cómo diferentes elementos culturales y patrimoniales se conectan en el espacio de forma que cualquiera lo pueda entender.
3.3 Comparación datos de fuente alternativa
Comparamos el número de museos (y monumentos) que hemos obtenido con la información obtenida de VisitValencia: “La cultura en València forma parte de la esencia del ciudadano manifestándose en más de 60 contenedores culturales entre museos, monumentos y espacios multidisciplinares.”
nrow(todos_museos)[1] 25La diferencia entre los 36 museos encontrados en nuestro conjunto de datos y los más de 60 mencionados por VisitValencia puede deberse a varias razones. En primer lugar, nuestro conjunto de datos puede estar incompleto o desactualizado, ya que solo se han incluido museos y monumentos específicos, mientras que VisitValencia contabiliza una mayor variedad de espacios culturales, centros multidisciplinares. Además, los criterios de inclusión pueden variar, ya que VisitValencia incluye tanto museos establecidos como exposiciones temporales, lo que podría aumentar el número total. También es posible que nuestra fuente de datos tenga una cobertura geográfica más limitada o parcial (sobre todo teniendo en cuenta que es de generación voluntaria) en comparación con la fuente oficial, lo que contribuye a la diferencia en los números.
Para obtener una representación más precisa de la cantidad total de museos y centros culturales en Valencia, sería recomendable ampliar y actualizar nuestra fuente de datos para incluir una gama más amplia de espacios y revisar los criterios utilizados en la recopilación de datos.
4. Bibliografía
Valencia, V. (2020, November 21). Monumentos y museos que no te puedes perder. Visit Valencia.
Gisadminbeers. (2019, September 4). Equivalencia de códigos EPSG - Gis&Beers.
Race. (n.d.). Zonas de bajas emisiones en España | RACE.
Coordenadas geográficas de España - Latitud y longitud. (n.d.). Geodatos.